Dielectric breakdown of Mott insulators in dynamical mean-field theory 
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Using nonequilibrium dynamical mean-field theory, we compute the time evolution of the current 
in a Mott insulator after a strong electric field is turned on. We observe the formation of a quasis- 
tationary state in which the current is almost time-independent although the system is constantly 
excited. At moderately strong fields this state is stable for quite long times. The stationary current 
exhibits a threshold behavior as a function of the field, in which the threshold increases with the 
Coulomb interaction and vanishes as the metal-insulator transition is approached. 



Nonequilibrium phase transitions and nonlinear trans- 
port are becoming central issues in the study of strongly 
correlated systems. One of the most basic phenom- 
ena is the dielectric breakdown (destruction of insulat- 
ing states due to strong electric fields) In Mott 
insulators, the electron motion is frozen as a result of 
strong repulsive interactions Q, and in equilibrium the 
doping of carriers into a Mott insulator leads to inter- 
esting quantum states such as high- Tc superconductivity 
in two-dimensional materials. Thus, it is natural to ask 
how nonequilibrium carriers behave when electrons in a 
Mott insulator start to move in response to strong electric 
fields. 

Experimentally, the physics of nonlinear transport in 
correlated electron systems has been studied in oxides P 
as well as in organic materials [Hj. In one dimensional 
Mott insulators, dielectric breakdown was observed and 
it was found that the current increases with a thresh- 
old behavior [ij. The current- voltage {I-V) characteris- 
tics exhibits a strong non-linearity with a negative dif- 
ferential resistivity between the weak current and large 
current regimes. More recently, the problem of nonlin- 
ear transport is attracting interest also in the cold atom 
community, where a novel realization of the Mott insu- 
lating state has been achieved. In Rcf. Q the effect of 
a potential gradient was studied to probe the excitation 
spectrum. 

A theoretical description of nonlinear transport is chal- 
lenging because one needs to take into account two non- 
perturbative effects, electric fields and electron-electron 
interactions, simultaneously. In one-dimensional sys- 
tems, reliable numerical techniques such as exact diago- 
nalization and the time-dependent density matrix rcnor- 
malization group are available, and a threshold behavior 
in the current was indeed observed Q]- A many-body 
Schwingcr-Landau-Zencr mechanism, in which doubly 
occupied states (doublons) and holes are pair produced 
by quantum tunneling was proposed as an explanation 
[3|. However, in previous studies, a direct calculation of 
the /-^-characteristics is absent since a steady state cur- 
rent cannot be easily reached in a finite system. This 
made it difficult to compare theory with experiments. 



Another totally unexplored issue is the temperature ef- 
fect, where experiments suggest a relatively strong tem- 
perature dependence of the current. 

The purpose of the present paper is to address these 
questions using dynamical mean-field theory (DMFT) 
[7|, which is suitable for the study of high-dimensional 
bulk systems. Nonequilibrium DMFT Q has been used 
to reveal various types of steady states and relaxation 
phenomena in the Falicov-Kimball model and 
the Hubbard model [l3| . For the current analysis one 
must compute the dynamics of the Hubbard model at 
rather strong interactions up to relatively long times. So 
far, this task has been prohibitively difficult for imp urity 
solvers based on real-time quantum Monte Carlo [13|, 



but it has become accessible recently through an imple- 
mentation of the self-consistent hybridization expansion 
within the Keldysh framework |14l |. 

In the following we focus on the Mott insulating phase 
in the half-filled Hubbard model on a d-dimensional cubic 
lattice with lattice spacing a. 



(1) 



Here {cia) denotes the creation (annihilation) oper- 
ator for an electron with spin a at lattice site Ri, U 
is the local Coulomb repulsion, and Vij describes hop- 
ping between the sites. To study the dielectric break- 
down of the Mott insulator, we initially prepare the sys- 
tem in thermal equilibrium at temperature T ~ 1/(3 
and apply a spacially homogeneous electric field F{t) 
for time t > 0. Using a gauge with pure vector po- 
tential A{t), i.e., F{t) ~ —dtA{t)/c, F{t) is incorpo- 
rated into Eq. ([T]) by means of the Peierls substitution, 
V,j{t) = cxp[ie(i2j - Ri)A{t)lhc]. The field is chosen 
to point along the body diagonal i) = (1 ... 1)* of the 
unit cell. It is turned on to a value F within a switch- 
ing time to, F[t) = riFr{t/to), using a switching profile 
r{x) = \- |cos(7ra;) -I- i cos(7rx)^ for < a; < 1. The 
main results turn out to be independent of the switching, 
and we choose io = 3 if not stated otherwise. 

In the limit d = oo the problem is solved exactly us- 
ing nonequilibrium DMFT. The DMFT self-consistency 



2 




2 4 6 8 10 12 14 0.5 

t F 



FIG. 1: Time-evolution after turning on a field F in the 
insulator {U = 5, /3 = 10): (a) The current j{t); results for 
F = 10 are multiplied with a factor 0.005 to match the scale, 
(b) Effective temperature Taft and (c), double occupancy d{t) 
for F = 0.8. The thin line shows the thermal expectation 
value of d at temperature T = Tcf!{E{t)). (d) Current at 
F = 0.8 for various switch-on times to- (e) Current, averaged 
for t > 10, obtained using either NCA or OCA to solve the 
DMFT impurity problem. 



equations for finite electric field are identical for the Hub- 
bard model and the Falicov-Kiniball model, and they are 
detailed in Ref. The DMFT single-site problem, 

on the other hand, is identical to the zero field case, 
and we solve it by means of the self-consistent hybridiza- 
tion expansion jl4| up to second order (one-crossing ap- 
proximation), which is known to be reliable for the in- 
sulating phase and the crossover regime Bin. We 
choose units of energy such that the density of states 
is given by p{e) oc exp(— e^), with a width (second mo- 
ment) W = V2/2. Time and field are measured in units 
of l/\/2W and \/2W/ea, respectively (h = 1). In equilib- 
rium, a first-order Mott transition is found in the para- 
magnetic phase 0, with a critical end-point at C/c ~ 3.1 
and Tc ~ 0.02 (determined within the one-crossing ap- 
proximation). 

After an electric field is turned on in the Mott- 
insulating phase, we either observe the formation of 
a quasistationary state with time-independent current 
(Fig.[T^), or, for very large values of F, the emergence of 
27r/i^-periodic Bloch oscillations {F = 10 in Fig. [T^). In 
analogy to the Falicov-Kimball model [13], where Bloch 
oscillations are quenched by the interaction, there is no 
sharp separation between oscillatory and non-oscillatory 
regimes, but instead the current behaves irregularly at in- 
termediate fields. In the following, we only study the qua- 
sistationary state, which will reveal the dielectric break- 
down of the Mott insulator at moderately large, possibly 
experimentally accessible fields. 



Since our system is not coupled to a thermal bath, the 
energy E = {H{t)) increases at a rate E = F{t)j{t), 
and thus a stationary state with nonzero current can- 
not exist forever. However, we find that j{t) remains 
remarkably stable even after a considerable energy in- 
crease. This fact becomes clearly evident if one looks 
at the effective temperature instead of the energy, i.e., 
the temperature Tf,ff(E) after a hypothetical thermaliza- 
tion at energy E: For F = 0.8 in Fig. [3 e.g., T^g in- 
creases by a factor 1.5 during times 4 < t < 15 (Fig.[Tl3), 
while j{t) remains almost constant. The saturation of the 
current and a simultaneous increase in the double occu- 
pancy d = {ni-fTiii) (Fig. [T]:) indicate that the current 
flow causes excitations in the system that are immobile, 
just like the spin fluctuations in the ground state that 
lead to a flnite double occupancy but no linear response 
conductivity. These excitations do not thcrmalize on the 
timescale of our simulation. Otherwise the current would 
increase thermally, and d{t) would have to match its ex- 
pectation value in thermal equilibrium at temperature 
Tcs{E{t)) (thin line in Fig. [T];), which is not the case. 
This behavior is consistent with recent experiments on 
ultracold gases [l^, where it was found that artificially 
created double occupancies in the Mott insulator relax 
only on exponentially long timescales. 

The quasistationary current turns out to be more or 
less independent on how the field is turned on, and by 
increasing the switching time to one can reduce the other- 
wise rather strong transient oscillations (Fig. [TJi). Even 
for slow switching, however, the transient current can be 
orders of magnitude larger than the stationary current. 
This fact is already entailed in the linear response rela- 
tion j(t) ~ J^ds a{s)F{t—s) which always holds for small 
enough times. The transient is at least proportional to 
F, while the long-time limit can be exponentially small 
(see below). In particular, in the Mott insulator at T = 
the dc conductivity ddc = /q dsa{s) vanishes, while the 

integral a = J^dt J^dsa{s) over the current yields a 
nonzero static polarizability. 

After averaging over times t > 10, the quasistationary 
current j at C/ = 5 shows a sharp increase around = 0.5 
(Fig. Wii), which is the hallmark of the dielectric break- 
down. We note that both first order (NCA) and second 
order (OCA) implementations of the self-consistent hy- 
bridization expansion yield similar results, but in anal- 
ogy to equilibrium calculations the insulating behavior 
is overestimated by NCA, such that the increase of the 
current is shifted to stronger fields. In the following we 
stick to the more reliable OCA as an impurity solver. 

A plot of the conductance j/F on a logarithmic scale 
reveals a crossover from the temperature-dependent lin- 
ear response conductivity (Tdc(T') at small F to an al- 
most temperature- independent curve at large F (Fig. [21a.). 
These data suggest that the quasistationary current has 
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FIG. 2: (a) Time- averaged current, ]/F, for f/ = 5: Cross- 
like and open symbols correspond to a time average for 8 < 
t < IQ and 10 < f < 12, respectively. Filled symbols result 
from an extrapolation to t = oo [see panel (b)]. The solid 
black line is obtained by fitting the data for /3 = 20 with 
Eq. [Fth = 1.92, (JiSn = 0.066]. (b) Long-time evolution 
of the current for f/ = 5, analyzed by a nonlinear least square 
fit j{t) — j + Sj exp(— t/r) (solid lines); the value J is shown 
by the filled symbols in (a), (c) The relaxation time r. 



a nonzero T — limit, which reflects the ground state 
decay, or dielectric breakdown of the insulator. In anal- 
ogy to the one-dimensional case 0, S], we will refer to 
this limiting value as the tunneling current j'tun- In fact, 
the low-temperature data in Fig. [2^ can be fit with the 
same exponential law that determines the ground state 
decay rate in the one-dimensional Hubbard model 0, Q 



Jtu„(F) = F<„exp(-Fth/F), 



(2) 



with a threshold Fth (solid line in Fig. Hi). 

For small fields and T > 0, the exponentially small 
Jtun(^) is dominated by the linear response current (JdcF. 
Surprisingly, however, we find that thermal and tunnel- 
ing current do not simply add up in the stationary state, 
but j/F can become much smaller than tTdc- The peculiar 
minimum at the crossover between tunneling and linear 
response regime in Fig. [2K arises because the relaxation 
to the stationary state becomes slower with decreasing 
field, such that an average of j(i) in a fixed time interval 
still yields adcF for — > 0. A detailed analysis of the 
long-time behavior at intermediate F reveals that the sta- 
tionary current is approached via an exponential decay 
j{t) ~ j + 5j exp(— t/r), where j is closer to jtun{F) than 
it is to the linear response current (Fig. [^b)- The decay 
of the thermal current indicates that mobile, thermally 




FIG. 3: (a) Linear response conductivity, obtained by extrap- 
olating j/F to F = 0. Thin solid lines correspond to fits with 
Eq. ([21; the resulting gap A{U) is shown in (b). (c) Conduc- 
tance j/F for the stationary current at /3 = 20. For U = 3.5, 
there is a still a slight drift of the current at the largest times 
(see text), and we plot time-averages for 8 < t < 10 (crosses) 
and 12 < f < 14 (open symbols). Solid lines are linear fits 
according to Eq. (d) The threshold field Fth(C/) resulting 
from the fits in (c). 



excited carriers (and those excited during the switch-on) 
arc transformed into immobile excitations in the presence 
of an electric field, such that the linear response current 
is only visible on a timcscalc cx 1/F. Clearly, such a be- 
havior must be specific to a closed system, whereas in an 
open system the coupling to the environment constantly 
tends to restore the thermal state. Nevertheless, the ob- 
served relaxation phenomenon is an interesting topic for 
future theoretical transport investigations, and it is an 
open question whether it is a generic feature in closed 
systems. 

Figure [3] summarizes the main numerical results of this 
paper by comparing the current in the tunneling regime 
(T I 0) and in the linear response regime — J- 0, t < 
1/-F). For t/ < 3, the linear response conductivity CTdc 
increases with decreasing temperature, while it becomes 
exponentially small for U >3, 



<7dc ^ exp(-A/T), 



(3) 



thus signaling the metal-insulator transition at /7 w 3 
(Fig. [5^ and b). (Strictly speaking, the metal- insulator 
transition displayed in Fig. [3] is only a crossover because 
temperatures are larger than Tc.) On the other hand, by 
fitting j/F for /3 = 20 with Eq. we obtain the thresh- 
old Fth as a function of the interaction (Fig. [3t and d). 
Deep in the insulating phase, both A and i^th increase 
linearly with U, while they vanish as the metal-insulator 
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transition is approached. Because the metal-insulator 
transition is first order in DMFT Q, observablcs in the 
insulating phase are expected to display nonanalytic be- 
havior not at the actual transition 11^.2 at T = 0, but 
at the spinodal point Uci where the insulator at T = 
disappears. However, the precise behavior of A and Fth 
at the transition is difficult to obtain within the current 
approach: As J^th decreases, a fit oi]/F with Eq. ([2]) be- 
comes increasingly difficult because the stationary state 
is eventually no longer reached within the numerically 
accessible times for F < Fth due to the slow saturation 
of j{t) for small F (cf. Fig. [It). 

Apart from the crossover region at U < 4, however, 
we do observe the behavior described by Eqs. ([2|) and ([3]) 
over a wide range of temperatures and fields. Although 
our setup involves a closed system in which energy is 
conserved, we expect that these DMFT results resemble 
the /-^-characteristics of a real Mott-insulating material 
for large and small fields, respectively. While the quasis- 
tationary state cannot exist forever in our setup, it can 
indeed be stable in a real solid provided the excess en- 
ergy is passed to the lattice. Since the tunneling current 
turns out to be essentially independent of the excitation 
of the system (cf. Fig. [TJ , differences between the quasis- 
tationary current in our setup and the stationary current 
in experiments are expected to be of the order of the lin- 
ear response current ([3]) , which is negligible compared to 
Eq. ([2]) for large enough fields. For small fields or high 
temperatures, on the other hand, the open system should 
recover the linear response behavior which is lost in the 
closed system at long times. Only in between the linear 
response and the tunneling regime one may thus expect 
that the current depends essentially on details that are 
not accounted for in the purely electronic description. 

We find that the ratio Fth /A is only weakly depen- 
dent on the interaction or the bandwidth in the insula- 
tor, and it should thus give the correct order of magni- 
tude for the breakdown field in real Mott insulators as 
well. Note that this value, i.e., Fth ~ 2 — 3 A/ea, is much 
larger than the temperature-dependent threshold which 
is obtained in the experiments of Ref. [l| in connection 
with a negative differential resistance (for SrCuOa, e.g., 
Fth ~ 10"'' A/ea at T = 190K). The threshold behav- 
ior in these one-dimensional materials must thus be of 
a different origin, and indeed collective excitations were 
proposed in Ref. as the temperature dependence of 
the experimental threshold is similar to what is expected 
for charge-ordered materials. The larger threshold found 
in our analysis should be observed in paramagnetic Mott 
insulators when other sources of destabilizing the insula- 
tor are not present. 

In conclusion, we have investigated the dielectric 
breakdown of a Mott insulator in the Hubbard model by 
computing the current in a strong electric field F. Our 
main result is the formation of a quasistationary noncqui- 
librium state with time- independent current, which may 



be called a field-induced metal. In the limit of small 
temperature, the stationary current resembles the expo- 
nential law [Eq. ([2])] for the ground state decay-rate in a 
one-dimensional Hubbard model due to many-body Lan- 
dau Zener tunneling [2], Q . Its value becomes exponen- 
tially small below a threshold field which vanishes at the 
metal-insulator transition. Although the threshold field 
is generally quite large it should be experimentally acces- 
sible, e.g., using thin films of insulating material between 
metallic leads. 
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